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In this paper, by using a stochastic reaction-diffusion-taxis model, we analyze the picophyto- 
plankton dynamics in the basin of the Mediterranean Sea, characterized by poorly mixed waters. 
The model includes intraspecific competition of picophytoplankton for light and nutrients. The 
multiplicative noise sources present in the model account for random fluctuations of environmental 
variables. Phytoplankton distributions obtained from the model show a good agreement with ex- 
perimental data sampled in two different sites of the Sicily Channel. The results could be extended 
to analyze data collected in different sites of the Mediterranean Sea and to devise predictive models 
for phytoplankton dynamics in oligotrophic waters. 
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I. INTRODUCTION 

Natural systems are characterized by two factors: (i) non-linear interactions among their parts; (ii) external pertur- 
bations, both deterministic and random, coming from the environmenlP^. It is worth noting that natural systems, 
because of these characteristics, are complex systems^ Therefore, the study of a marine ecosystem has to be 
performed by considering the perturbations, not only deterministic but also random, due to the fluctuations of the 
environmental variables. This implies the necessity of including in the model a term which describes the continuous 
interaction between the ecosystem and environment. In particular, physical variables, such as temperature, salinity 
and velocity field, are affected by random perturbations and can be therefore treated as noise sources. This causes 
the phytoplankton behaviour to be subject to a stochastic dynamics, and allows to expect that a stochastic approach 
should reproduce the distributions of phytoplankton biomass better than deterministic models. On this basis, noise 
effects have to be included to better analyze the dynamics of a marine system such as that studie d in this work. 
The growth of phytoplankton is limited by the concentration of nutrients R and intensity of light In particular, 

the survivance of phytoplankton is strictly connected with the presence of sufficiently high nutrient concentration. 
It is worth stressing that nutrients, which are in solution, diffuse from the bottom (seabed) towards the top (water 
surface). Nutrient distributions along the water column are therefore characterized by an increasing trend from the 
sea surface to the benthic layer. As a consequence, the positive gradient of nutrient concentration causes the maxima 
of chlorophyll, which is contained in the phytoplankton cells, to be localized in deep subsurface layers. This condition 
constitutes one of the most striking feature of the nutrient poor waters in ocean ecosystems and freshwater lakeS-^-"^. 
Conversely, the light penetrates through the surface of the water and has an exponentially decreasing trend along 
the water column. This characteristic makes the deep layers unfavourable for the photosynthesis, determining, as a 
consequence, adverse life conditions for phytoplankton. In particular, light is a crucial parameter for the localization 
of the deep chlorophyll maximum (DCM), as revealed by the significant correlation found between the depth of DCM 
and light intensity over the Mediterranean basin in summer (Brunet et al., unpublished data). 

The dynamics, competition and str ucturing of phytoplankton populations have been investigated in a series of theoret- 
ical studies based on model systemPSEIEMSIl, jn a few recent investigations it was observed that in the presence of an 
upper mixed layer either surface or deep maxima can be observed indifferently under almost the same conditions'^^ 
In view of analyzing an ecological system, as a preliminary step it is necessary to define the correct values of the pa- 
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rameters and the role that they play on the dynamics of the populations, specifically when the coexistence of different 
species in the same community is considered^^. The responses of the species to environmental solicitations strongly 
depend on the biological and physical parameters. Among these, a relevant role is played by the phytoplankton 
velocity which is strictly connected with the microorganism size, one of the main functional traits for phytoplankton 
diversity. Oth er par ameters that influence the balance of a marine ecosystem are, for example, growth rates and 
nutrient uptak^Sfilsg] 

In this paper we deal with data obtained in a hydrologically stable area of the Mediterranean Sea, where the environ- 
mental light and nutrients, specifically phosphorus, contribute to determine life conditions. The Mediterranean basin 
is characterized by oligotrophic conditions and it has been suggested that there is a decreasing trend over time in 
chlorophyll concentration. This has been associated with increased nutrient limitation resulting from reduced vertical 
mixing due to a more stable stratification of the basin, in line with the general warming of the Mediterranean 
Here we consider the Strait of Sicily, which is known to govern the exchanges between the eastern and western basins 
and is characterized by active mesoscale dynamics^^, strongly influencing the ecology of phytoplankton communities. 
Mor eover , the Strait of Sicily is a biologically rich area of the Mediterranean Sea with a key role in terms of fish- 
erie^^^IlSl^ The anchovy growth (along with phytoplankton biomass) in the Sicilian Channel resulted to be mainly 
explained by changes in the chlorophyll concentration, used as a phytoplankton biomass indicator Our study is 
performed using a stochastic model obtained by modifying a deterministic reaction-diffusion-taxis model. Specifically, 
the analysis focuses on the spatio-temporal dynamics of the phytoplankton biomass, and provides the time evolution 
of biomass concentration along the water column. Finally, the results are compared with experimental data collected 
in two different sites of the Strait of Sicily. 

II. MATERIALS AND METHODS 
A. Environmental data 

The experimental data were collected in the period 12*'' - 24*'* August 2006 in the Sicily Channel area (Fig. [ij 
during the MedSudMed-06 Oceanographic Survey onboard the R/V Urania. Hydrological data were obtained using a 
SBE911 plus CTD probe (Sea-Bird Inc.); chlorophyll a fluorescence data {chl a, /ig/1) were contemporary acquired by 
means of the Chelsea Aqua 3 sensor. In the Libyan area the CTD stations were located on a grid of 12 x 12 nautical 
miles. Moreover, CTD data have been collected along a transect between the Sicilian and the Libyan coasts. In the 
present work, two stations out of the whole data set were considered. The selected stations were located on the south 
of Malta (site L1105) and on the Libyan continental shelf (site L1129b). The collected data were quality-checked and 




FIG. 1: Locations of the CTD stations where the experimental data were collected. 
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processed following the MODB instruction^^ using Seasoft software. The post-processing procedure generated a text 
file for each station where the values of the oceanographic parameters were estimated with aim step. Hydrological 
conditions remained constant for the entire sampling period and were representative of the oligotrophic Mediterranean 
Sea in summer. Nitrate, nitrite, silicate and phosphate concentrations were not determined. 



B. Phytoplanktonic data 

Depending on size the phytoplankton species can be divide into two main fraction: 

• < 3/im picophytoplankton, formed by groups, Prochlorococcus, Synechococcus and picoeukaryote^^^'^. This 
size of phytoplankton accounts for about 80% of the total chl a on average'^, ranging from 40% to 90% (69% 
in the DCM)49. 

• > 3/zm nano- and micro-phytoplankton, characterized by a lower correlation with nutrients and salinity respect 
to picophytoplankton. This is connected with the fact that the contribution of picophytoplankton in the DCM 
is higher than in the surface layei'^. This larger size fraction of phytoplankton amounts to 20% of the total chl 
a on average and is uniformly distributed along the water column. 

The high pigment diversity of the smaller phytoplankton in the DCM and its elevated contrib ution to the total chl 
a indicated a strong degree of adaptation to the quantity and quality of light availablj ^^ l ^^ l ^^ l This is not true for 
the larger phytoplankton, which is represented mainly by diatoms or Haptophytes. Picoeukaryotes, which belong 
to the smaller size class, present peculiar eco-physiological properties-*^----^'^ , such as low sinking, high growth rate 
and low nutrient uptake. Their small size leads to a low package eff ect, w hich contributes to the light-saturated rate 
of photosynthesis that can be achieved at relatively low irradiance^^^ElHIll, to their peculiarities and relevant 

role in ecosyste m functio ning, they constitute a key-group to be considered within a model of population dynamics. 
In Sicily ChanneP^'^SllII^ picophytoplankton is numerically dominated by the Prochlorococcus fraction. In this area 
the number of Prochlorococcus cells is constant in the first 20 m, and is characterized in the DCM by an average 
value of 5.2 x 10^ cell ml~^. Average picoeukaryote concentration in the DCM is 0.6 ± 0.4 x 10^ cell ml~^, and the 
mean value of chl a cell~^ ranges between 10 and 660 fg chl a cell""'^ along the water column, with a significant 
exponential increase with depth (see Fig. [2p2!. The concentration of chl a (fg cell"^) per cell in picoeukaryotes was 
highly variable among different water masses, with signifi cantly hig her values in the DCM respect to the surface, as 
a result of photoacclimation to decreased light irradiance d'^^l^" * ^^ -^. 
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FIG. 2: Mean vertical profile of chl a per picoeukaryote cell (fg cell~^). Error bars are Standard Deviation. Equation and 
for the fit are reported on the plots. (Courtesy of Brunet et al., 2007). 
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III. EXPERIMENTAL RESULTS 

Data obtained from the cruises in two different sites of the Strait of Sicily both for temperature and chl a con- 
centration are shown in Fig. [3j In site L1129b, the behaviour of the temperature along the water column indicates 
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FIG. 3: Profiles of temperature (panels a, c) and chl a concentration (panels b, d) measured in sites L1129b and L1105. The 
black lines have been obtained by connecting the experimental points corresponding to samples distanced of 1 meter along the 
water column. The total number of samples measured in the two sites is n = 176 for L1129b, and n = 563 for L1105. 

the presence of a mixed layer (from the surface to 28 m depth) characterized by a high value of temperature. Below 
the thermocline (28 m depth) the temperature decreases up to 80 m, becoming uniform below this depth (Fig. [3^). 
The site L1105 shows a mixed layer over the first 24 m of depth, and a sharp decrease of temperature from 24 to 75 
m (Fig. [sj;). Experimental data for chl a concentration show a nonmonotonic behaviour, as a function of the depth, 
characterized by the presence of DCM in both sites (see Fig. [3]3,d). Specifically, fluorescence profiles show a similar 
behaviour in the two sites, with chl a concentration ranging between 0.010 and 0.17^g chl a 1"^. Differences between 
the two sites are observed in the depth, shape and width of the DCM. 

IV. THE MODEL 

In this study we analyze the spatio-temporal dynamics of a picophytoplankton community, limited by nutrient 
and light in a vertical poorly mixed water column. The mechanism, responsible for the phytoplankton dynamics, is 
schematically shown in Fig. 4. The mathematical tool used to simulate the phytoplankton dynamics is an advection- 
reaction-diffusion model. In particular, we investigate the distribution of the picophytoplankton along the water 
column, with light intensity decreasing and nutrient concentration increasing with depth. Analysis and numerical 
elaborations are divided in two phases: 

• Phase 1. By using a model based on two differential equations, the distribution of picophytoplankton biomass 
h is obtained along the poorly mixed water column as a function of the time and depth, and simultaneously the 
distribution of nutrient concentration i?, which limits the growth of phytoplankton, is calculated. The results 
obtained are compared with the experimental data collected in the two different sites of the Strait of Sicily. 
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FIG. 4: (Color online) Scheme of the mechanism responsible for the phytoplankton dynamics (modified from original figure by 
A. Ryabov). (a) Image of Micromonas N0UM17 (courtesy of Augustin Engman, Rory Welsh, and Alexandra Worden. 



• Phase 2. In order to match better the results for b and R to the experimental data, the random fluctuations 
of the environmental variables are taken into account. In particular, a stochastic model is obtained from the 
deterministic one by inserting into the equations terms of multiplicative Gaussian noise. 



A. The deterministic model 



Here we introduce the model consisting of a system of differential equations, with partial derivatives in time 
and space (depth). The model allows to obtain the dynamics of the phytoplankton biomass b{z,t) and nutrient 
concentration R{z,t). The light intensity I{z,t) is given by a function varying, along the water column, with the 
depth and biomass concentration. The behaviour of the phytoplankton biomass, along the water col umn, is the re sults 
of three processes: growth, loss, and movement. The phytoplankton growth rate depends on / and J j20 | 2i | 32 | 33 | 36 | r^^ie 
limitation in phytoplankton growth is described by the Monod kinetic^^. The gross phytoplankton growth rate per 
capita is given by min{//(/), //j(i?)}, where //(/) and ,fri{R) are obtained by the Michaelis-Menten formulas 

fiiI)^rI/iI + Kj), (I) 
MR)^rR/iR + Kn). (2) 

In Eqs. (1), (2), r is the maximum growth rate, while Kj and Kn are the half-saturation constants for light intensity 
and nutrient concentration, respectively. Varying Kj^ and Kj allows to model, for instance, a species which is better 
adapted to the light (smaller values of Kj) or nutrient (smaller values of Kfj). More specifically, we consider a species 
with small Kj and large Kr that corresponds to good life conditions at large depth. These constants depend on the 
metabolism of the specific microorganism considered. 

The biomass loss, connected with respiration, death, and grazing, occurs at a rate jtJMMHI, xhe gross per capita 
growth rate is defined as 



giz,t) ^ mm{fR{R{z,t))Jj{I{z,t))). (3) 

Turbulence, responsible for passive movement of the phytoplankton, is modeled by eddy diffusion. Specifically, 
we describe turbulence assuming that the vertical diffusion coefficient is uniform with the depth and characterized 
by a low value {Di, = Dr = 0.5). This choice is motivated by the fact that in sites L1129b and L1105 the 
phytoplankton peaks, located at 87 m and III m respectively, are quite far from the thermocline (see Fig. |3|. 
Therefore, phytoplankton should go up (or down) if the biological conditions are more suitable for growth above 
(below) than below (above). Finally, no migration should occur if the biomass concentrations are the same at 
different depths. These assumptions about growth, loss, and movement, allow to obtain the following differential 
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equation for the dynamics of biomass concentration pHMl 



db{z,t) 

dt 



^ d'^b(z,t) db(z,t) 
g{z,t)b{z,t) ~ mb{z,t) + Db „ ^ ' -v- ^ ' ' 



dz 



dz 



(4) 



The positive phytoplankton velocity v, due to active movement, is oriented downward (sinking), in the direction of 
positive z. Phytoplankton does not enter or leave the water column. This is set by using no-flux boundary conditions 
at z = and z = Zf,: 



db 

dz 



db 
dz 



= 0. 



(5) 



Eddy diffusion is responsible for mixing of the nutrient concentration along the water column, with diffusion co- 
efficient Dfj. The nutrient consumed by the phytoplankton is also obtained from recycled dead phytoplanktonic 
microorganisms. The dynamics of nutrient concentration can be therefore modeled as follows: 



dR{z,t) 
dt 



h{z,t) 
Y 



giz,t) + DR 



d^Rjz^t) 
dz'^ 



em- 



h {z,t) 
Y '■ 



(6) 



Here Y is the phytoplankton produced biomass per unit of consumed nutrient, and e is the nutrient recycle coefficient. 
Since the nutrient is not supplied by the sea surface but comes from the seabed, its concentration is set to the 
constant value Rin in the sediment and, as a consequence, to the value i?(z&) in the bottom of the water column. In 
fact the nutrient diffuses across the sediment-water interface with a rate proportional to the concentration difference 
between the solid phase (seabed) and the deepest water layer (bottom of the water column). Accordingly, the 
boundary conditions are given by: 



= h{R,,,- R{zb)), (7) 

where h is the permeability of the interface. Finally, taking into account Lamber-Beer's lawE^Ml^ ^j^g light intensity 
is characterized by an exponential decrease modeled as follows: 

J(z) = /„exp|-^ [ab{Z) + abg\dZ^, (8) 

where a and and a^g are phytoplankton biomass and background attenuation coefficients, respectively. Equations 
Q-Q form the biophysical model used in our study. 



dR 

dz 



0, 



dR 

dz 



B. Results of the deterministic model 



The time evolution of the system is studied by analyzing the spatio-temporal dynamics of biomass and nutrient 
concentrations. In particular, by using a numerical method, implemented by a program in CH — h language and based 
on an explicit finite difference scheme, equations Q-Q are solved. The increment of the spatial variable is set to 
0.5 m. In view of reproducing the spatial distributions observed in the real data for the phytoplankton biomass (see 
Fig. |3| , we choose the values of the environmental and biol ogica,l p arameters to satisfy the monostability condition 
corresponding to the presence of a deep chlorophyll maximunPSESiSai^ Xhe numerical values assigned to the parameters 
are shown in Table |l] Specifically, the values of the biological parameters r, Kj, K^, v, have been chosen to reproduce 
the behaviour of picoeukaryotes. We note that, in systems characterized by a constant value o f the diffusion coefficient, 
the stationary state does not depend on the initial conditions, according to previous studieJ^SEH. in order to obtain 
the steady spatial distribution, we integrated numerically our equations over a time interval long enough to observe 
the stationary solution. As initial conditions we consider that the phytoplankton biomass is concentrated in the 
layer where the maximum of the experimental chlorophyll distribution is observed. On the other side the nutrient 
concentration is approximately constant from the water surface to the DCM, and increases linearly below this point 
up to the seabed. 

Preliminary analysis (data not shown) revealed that the stationary solution is characterized by DCMs which are 
shallower as the nutrient supply increases, and deeper for enhanced light radiation. In general, large values of 
lin (incident light intensity at the water surface) lead to stationary conditions characterized by DCM, while large 
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Symbol 


Interpretation 


Units 


Site L1129b 


Site L1105 




Incident light intensity 


/iuiol photons m~^ 


1404.44 


1383.19 


abg 


Background turbidity 




0.045 


0.045 


a 


Absorption coefficient of phytoplankton 




6 X 10"'° 


6 X 10""' 


Zb 


Depth of the water column 


in 


186 


575 


Db = Dr 


Vertical tmbulent diffusivity 




0.5 


0.5 


r 


Maximum sjiiecific growth rate 


h-1 


0.08 


0.08 


K, 


Half-saturation constant of light-limited growth 


/.miol jiihotons m~^ s"-'- 


20 


20 


Kb 


Half-saturation constant of rnitr lent- limited growth 


mmol nutrient m~'^ 


0.0425 


0.0425 


m 


Specific loss rate 


h-i 


0.01 


0.01 


IIY 


Nutrient content of phytoplankton 


mmol nutrient cell"-^ 


1 X 10"^ 


1 X 10"" 


€ 


Nutrient recycling coefhcient 


dimensionless 


0.5 


0.5 




Buoyancy velocity 




—0 0042 


—0 0042 




Nutrient concentration at Zf) 


mmol nutrient m~'^ 


26.0 


36.0 


h 


Sediment-water colunm permeability 




0.01 


0.01 



TABLE I: Parameters used in the model. The values of the biological and environmental parameters are those typical of 
picophytoplankton and summer period in Mediterranean Sea, respectively. 

values of Rin (nutrient concentration in the sediment) determine an upper chlorophyll maximum (UCM). Finally, 
for intermediate values of lin and Rin the chlorophyll maximum can be localized close to the surface or at different 
depths, depending on the values of the other parameter^SH^ 

In our study the values of the light intensity resulted to be quite high in both sites, since sampling occured during 
summer (August 2006). In this period the light intensity at the water surface is larger than 1300 /zmol photons 
m~^ s~^. Moreover the sinking velocity is set to the value typical for picophytoplankton, v ~ 0.1 m day"-''^. The 
diffusion coefficent is fixed at the value Di, = 0.5 cm^ sec~^, which corresponds to the condition of poorly mixed 
waters. By solving Eqs. we obtain the biomass concentration expressed in cells/m'^ along the water column. 

Depths of the water column used in the model were set according to the measured depths in the corresponding marine 
sites. Moreover the light intensities, are fixed using data available on the NASA web sit^Sl!_ Finally, nutrient 
concentrations at the seabed were set at values such as to obtain, for each site, a peak of biomass concentration at the 
same position of the peak experimentally observed. All the other parameters are the same in both sites. The growth 
rate obtained from Eq. ^ agrees with the values measured by other authord^^. 

We note that our numerical results were obtained using a maximum simulation time tmax = 10^ h. Simulations (here 
not reported) performed within the deterministic approach show that the stationary regime is reached at i « 3 • 10^ h. 
This indicates that, to reach the steady state, it is sufficient to solve the equations of our model with a maximum 
time t„iax = 4 • 10^ h. By this way, we get the stationary profiles, both for biomass concentration and light intensity, 
shown in Fig. [Sj Here we can note the presence of a biomass peak as found in the experimental data, and the typical 
exponential behaviour of the light intensity. To compare the theoretical results with the experimental data, we 
exploit the curve of Fig. [2] to convert the cell concentrations, obtained from the model and expressed in cell/m'^, into 
chl a concentrations expressed in fig/l. We recall that about 43% of the total quantity of chl a^s 59 is due to nano- 
and micro-phytoplankton (20% of the total chl a on average), and Synechococcus (23% of the total chl a on average), 
quite uniformly distributed along the water column. Since our model accounts for the dynamics of picoeukaryotes, to 
compare the numerical results with the experimental data, we consider the 43% of the total biomass and divide it by 
depth, obtaining for each site the value Abdu a , which represents a constant concentration due to other phytoplankton 
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FIG. 5: Stationary distributions of the biomass concentration and light intensity in sites L1129b (panels a, b) and L1105 (panels 
c, d) as a function of depth. 




FIG. 6: (Color online). Stationary distribution of the chlorophyll a concentration as a function of depth calculated (red line) 
by the deterministic model and measured (green line) in sites (a) L1129b and (b) L1105. 



species present in the water column. Finally along the water column we add the theoretical concentration with /Sbcu a 
and obtain, for the distributions of chl a concentration, the stationary theoretical profiles consistent with those of 
the experimental data. The results are shown in Fig. |6] Here we can observe that in both sites the deep chlorophyll 
maxima obtained from the model are located at the same depth of those observed experimentally. However, the shape 
of the theoretical chl a distributions is quite different from the experimental profiles. Finally, we note that in site 
L1105 the magnitude of the theoretical DCM is significantly different from that observed in real data. 
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C. The Stochastic Model 



In the previous section we used a deterministic model to fit the experimental distributions of chl a concentration. 
The results obtained reproduce partially the characteristics of the experimental profiles. In order to get a good 
agreement between real data and theoretical results, we recall that the sea is a complex system. This implies, as 
discussed in Introduction, the presence of non-linear interactions among its parta^^Sl and a continuous interaction 
between the ecosystem and environment. In particular, the system dynamics is affected not only by deterministic 
forces but also random perturbations coming from the environment. In this context environmental variables, due to 
their random fluctuations, can act as noise sources, causing phytoplankton to be subject to a stochastic dynamics. 
Therefore, in order to perform an analysis that takes account for real conditions of the ecosystem, it is necessary to 
modify our model, including the noise effects. In the following we analyze two different situations. 
Case 1. The environmental noise affects only the biomass concentration. Therefore, Eqs. ([5])-([8]) are maintained 
unaltered, while Eq. Q becomes 



db , , ^ d^b db , ^ , , 

^=9b-mb + Db^-v—+b£,b{z,t) 9 
ot oz'^ oz 

Case 2. The environmental noise affects only the nutrient concentration. In this case, Eqs. Q , ([s]) , ([7| , (|8| are 
maintained unaltered, while Eq. (|6]) is replaced by 



dR 

Ik 



\me - 



g] 



Y 



'dz^ 



Riniz.t). 



(10) 



In Eqs. ([9]) and (10), ^fc(z,t) and ^ii[z,t) are statically independent white Gaussian noises with the usual properties 
(6(z,i)) = 0, = 0, {Uz,t)Uz',t')) = atSiz - z')5{t - t'), {U{z,t)iR{z' ,t')) = grS^z - z')5{t - t% where 

tTft and CTfl, are the noise intensities. We note that the two noise sources are spatially uncorrelated, that is at the 
generic point z no effects is present due to random fluctuations occurring in z' ^ z. 



D. Results of the stochastic model 



In this paragraph we solve numerically, within the Ito scheme, the equations of the stochastic model for different 
values of the noise intensities, obtaining the distributions of the picophytoplankton concentration as an average over 
1000 realizations. We recall that the ecosystem is characterized by non-linear interactions among its parts. Because 
of this feature the response of the system to external solicitations is also non-linear. Therefore, one can not expect 
that the presence of a symmetric noise with zero mean, i.e. Gaussian noise used in the model, produces in average 
the same effect as a deterministic dynamic^^. On the other side, the use of a random function, i.e. noise source, to 
simulate the spatio-temporal behaviour of the system, makes the single realization unpredictable and unique, and 
therefore non-representative of the real dynamics. As a consequence, one possible choice to describe correctly the 
time evolution of the system is to calculate the average of several realizations. This procedure, indeed, allows to 
take into account different "trajectories" obtained by the integrat ion of the stochastic equations, without focusing 
on a specific realization^ . According to the discussion of Paragraph IV B we calculated the solutions for a maximum 
simulation time tmax = 4 • lO** h. In Figs. [7]and|8]we show the results for case 1. Here we note that, in both sites, 
for higher noise intensities the peaks of the two average chl a distributions show: (i) a decrease of their magnitude; 
(ii) a small displacement along the water column. For suitable values of the noise intensity the peaks of the average 
chl a distributions obtained from the model match very well the experimental data. We note also that the two DCMs 
are located at 90 m (site L1129b) and 106 m (site L1105) of depth (in Figs. [7[i and |8|l compare theoretical (red 
line) and experimental (green line) profiles). To better understand the dependence of the biomass concentration 
on the random fluctuations of the nutrient, according to the procedure followed for case 1, we study for both sites 
the behaviour of the depth, width, and magnitude of the DCM as a function of gr. A quantitative comparison of 
each theoretical chl a distribution (red line) with the corresponding experimental one (green line) was carried out by 
performing goodness-of-fit test. The results are shown in Tables |llj where indicates the reduced chi-square. 
Results of the test show that the smallest difference between theoretical and experimental chl a distributions is 
obtained for ai, = 0.22 in site L1129b and 0-5 = 0.15 in site L1105. We also note that the depths of the DCMs are 
almost the same as in the deterministic case. 

In order to better analyze this aspect, we study for both sites the behaviour of the magnitude, depth, and width of 
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FIG. 7: (Color online). Average c/i/ a concentration calculated (red line) for different values of at by the stochastic model (case 
1, see Eqs. ([sjl, (|6|, ([7|, (|8|, ([ojl) as a function of depth. Results are compared with chl a distributions measured (green line) in 
site LI 129b. The theoretical values were obtained averaging over 1000 numerical realizations. The values of the parameters are 
those shown in Tablejl] The noise intensities are: (a) ai, — (deterministic case), (b) at — 0.10, (c) CT(, = 0.20, (d) ai, = 0.22, 
(e) at = 0.25 and (f) at = 0.30. 



the DCM as a function of crfc. The results, shown in Fig. |9j indicate that the depth of the DCM is almost constant 
for at < 0.4, increasing for higher values of the noise intensity (see panels b, e of Fig. |9]). Conversely, the width of 
DCM is characterized by a non-monotonic behaviour for increasing noise intensities. In particular, we note that the 
width of the DCM exhibits a maximum in both sites (for ah < 0.4 in site L1129b and at < 0.3 in site L1105). For 
higher noise intensities the width tends to zero for site L1129b, while a minimum is present for site L1105 at df, < 0.5. 
However, for at > 0.4, the values of the DCM width are less significant, since the chl a concentration along the water 
column and in particular in the DCM decrease strongly, as can be checked in panels a, d. In particular, random 
fluctuations, cause the reduction of biomass concentration and its displacement along the water column, determining 
the extinction of the picophytoplankton in the presence of higher intensities of noise. In this condition a clear 
determination of the DCM becomes more difficult. As a consequence, the values of depth and width for the DCM are 
less reliable. This analysis shows that the stationary conditions of the system depends strongly on the environmental 
fluctuations, which play a critical role in determining the best life conditions for the picophytoplankton species. 

We complete the analysis of the stochastic dynamics, considering the noise source which affects directly the nutrient 
concentration (case 2). By numerically solving the corresponding equations of motion (see Eqs. (|4|),([5|,([7l),([8|),(10)) 
and averaging over 1000 realizations, we obtain the average chl a distributions shown in Figs. [T0] and |ll| The results 
show that also for low noise intensities {an between 0.001 and 0.005), a decrease and a deeper localization of the 
DCMs are present. The shape of the chl a peaks exhibits, for both sites, a better agreement with the corresponding 
experimental DCMs respect to the deterministic case. In particular, for site L1129b the best value of the t^st 
is obtained for aj^ = 0.0020, while for site L1105 the best fitting results for aji = 0.0015 (see Table III). We note 
that in site L1129b the best agreement between experimental and numerical distributions is obtained, both in case 
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FIG. 8: (Color ordine). Average chl a concentration calculated (red line) for different values of ab by the stochastic model (case 
1, see Eqs. ([sjl, (|6|, ([7|, (|8|, ^) as a function of depth. Results are compared with chl a distributions measured (green line) in 
site L1105. The theoretical values were obtained averaging over 1000 numerical realizations. The values of the parameters are 
those shown in Table|l] The noise intensities are: (a) at — (deterministic case), (b) at = 0.05, (c) at = 0.10, (d) at = 0.15, 
(e) at = 0.20 and (f) at = 0.30. 

1 and case 2, for values of the noise intensity, ab and au, higher than those of site L1105. This can be explained by 
the fact that in site L1129b the DCM is localized at a depth shallower than in site L1105 (88m vs. Ill m), causing 
the environmental variables to be subject to more intense random fluctuations due to the closer sea surface. As 
a consequence, the chl a peak in site L1129b (88 m) is more strongly affected by the environmental noise than in 



site L1105 (111 m). The results, shown in Fig. 12 indicate that the depth of the DCM slightly increases in both 



sites as a function of the noise intensity (see panels b, e of Fig. 12). We note also that a decrease of the chl a 
concentration is observed in the DCMs of the two sites. This decrease is more rapid in site LI 105 (panel d), where a 
chl a concentration ^ 0.025 is reached for ~ 0.01. Analogously we observe an increase, faster in site L1105, of the 
width of the DCM. The spread of DCM and reduction of its magnitude are strictly connected with each other. In 
fact, the decrease of chl a concentration determines a flattening of the DCM with a consequent increase of its width. 
In conclusion the results shown in Fig. [12] indicate that the phytoplankton biomass tends to disappear for aji ~ 0.01, 
a value lower than those used in case 1, where no extinction occurs up to ab ~ 0.7 (see panels a, d of Fig. [9]). This 
indicates that the stability of the nutrient concentration is a critical factor in the dynamics of the ecosystem. Indeed, 
random fluctuations of the nutrient concentration can produce dramatic effects such as the collapse of phytoplankton 
biomass considered in our model. 

The previous analysis indicates that our model is able to reproduce the phytoplankton distributions observed 
in real data, without the model taking into account explicitly the environmental variables such as salinity and 
temperature. However, we observe that, in case 2, the spatio-temporal dynamics of nutrients has been modeled by 
introducing noise sources which can be interpreted as the effect of random fluctuations of environmental variables, 
among which salinity and temperature. 
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Site L1129b 



Site L1105 
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-r/2 
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T) 
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26 


0.00 


4.43 


0.0253 




36 


0.00 


22.87 


0.0407 


26 


0.10 


3.79 


0.0216 




36 


0.05 


22.72 


0.0404 


26 


0.20 


3.45 


0.0197 
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0.10 


22.55 


0.0401 


26 


0.22 


3.44 


0.0196 
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0.15 


22.50 


0.0400 


26 


0.25 


3.46 


0.0198 




36 


0.20 


22.95 


0.0408 


26 


0.30 


3.60 


0.0206 




36 


0.30 


27.14 


0.0483 



TABLE II: Results of x^, reduced chi -square (x^) goodness-of-fit test for site L1129b (left panel) and site L1105 (right panel) 
for different values of ai, (stochastic dynamics - case 1). The number of samples along the water column is n = 176 for site 



L1129b and n = 563 for site L1105. 
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Depth of DCM 



Width of DCM 





FIG. 9: Magnitude, depth, and width of the DCM as a function of (Jb obtained from the model for site L1129b (panels a, b, c) 
and site L1105 (panels d, e, f). 



V. DISCUSSION AND CONCLUSIONS 

In this work we presented a stochastic model, devised starting from previous deterministic model^^^Ml^ to study the 
spatio-temporal dynamics of the phytoplankton biomass along water column in two different sites of Sicily Channel. 
In our study, for fixed v, we chose values of the vertical turbulent diffusivity Dh which determine the absence of 
intrinsic oscillations of the phytoplankton concentration, maintaining the system far from the chaos. In oligotrophic 
waters, typical of Mediterranean Sea, where the surface mixed layer is depleted of nutrients, subsurface maxima of 
chlorophyll concentration and phytoplankton biomass are often found. Such deep chlorophyll maxima are permanent 
features in large parts of the tropical and subtropical oceans^^-S^-^. Furthermore, seasonal DCMs commonly develop 
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FIG. 10: (Color online). Average chl a concentration calculated (red line) for different values of an by the stochastic model 
(case 2, see Eqs. Q, Q, ([Tojl) as a function of depth. Results are compared with chl a distributions measured 

(green line) in site L1129b. The theoretical values were obtained averaging over 1000 numerical realizations. The values of 
the parameters are those shown in Table |l] The noise intensities are: (a) an — (deterministic case), (b) an = 0.0010, (c) 
cTfl = 0.0015, (d) o-jj = 0.0020, (e) an = 0.0025 and (f) an = 0.0050. 



in temperate regionP^'SZI and even in the polar ocean^^H^ when nutrients are dep leted in the surface layer with the 
onset of the summer season. Here we extend recent phytoplankton model^^SIinillHZIl to show that the phytoplankton 
distributions, due to random changes, can exhibit fluctuations. 

Our work consists in the analysis and subsequent modelhng, based on stochastic equations, of data from Sicily Chan- 
nel, where the waters are prevalently oligotrophic, the climatic conditions are those typical of a temperate region, and 
the DCMs show stable features for given conditions of hght and food resources. For values of depth ranging from 60 
to 110 meters the presence of a deep chlorophyll maximum indicates the existence of favourable life conditions for the 
phytoplankton and results in a good agreement with other experimental works, where higher biomass concentration 
and greater diversity are observed between 60 and 90 meters. At the depths considered in this work the light intensity 
is strongly reduced respect to the surface value (1% of the surface irradiance at 75 m). However, the low light intensity 
did not appear to limit the diversification of the phytoplankton communitj33sII_ jj^ fact, at depths ranging from 60 
to 90 meters a greater bio-diversity is observed. This can be explained considering that, at these values of depth, the 
high concentration of nutrients determines the most favourable life conditions for many species of phytoplanktoil^. 
Differences in the composition of phytoplankton between the surface and the DCMs are evident mainly for the smaller 
size class (less than 3 fim), which exhibits greater bio-diversity at depths between 60 and 90 meters. This could be 
due to the fact that different species of phytoplankton exhibit different responses to the limiting conditions. We recall 
that in the marine sites analyzed in this work the incident light intensity is characterized by high values (/„ > 1300 
fimol photon m~^s~^). Therefore, close to the surface the low nutrient concentration represents a limiting condition 
for all the phytoplankton species, so that the biomass concentration increases with depth. However, for larger values 
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FIG. 11: (Color online). Average chl a concentration calculated (red line) for different values of gr by the stochastic model 
(case 2, see Eqs. Q, Q, ([lOjl) as a function of depth. Results are compared with chl a distributions measured 

(green line) in site LI 105. The theoretical values were obtained averaging over 1000 numerical realizations. The values of 
the parameters are those shown in Table |l] The noise intensities are: (a) an ~ (deterministic case), (b) an = 0.0010, (c) 
cTfl = 0.0015, (d) o-jj = 0.0020, (e) an = 0.0025 and (f) an = 0.0050. 



of depth the light intensity becomes a main limiting factor for so me species, such as Synechococcus, which show a 
low degree of adaptability to smaller values of light intensitjE^. This causes Pr ochloroc occus and picoeukaryotes, 
which show a high degree of genetic plasticity ^ and tolerate lower light intensitie^^^EIlIIl^ exhibit a dominance in 
the deep chlorophyll maximum*^. 

In our model, the values of the biological parameters are those of the picoeukaryotes and the environmental pa- 
rameters are set at values typical of the oligotrophic waters during the warm period. These values allow to obtain chl 
a distributions along the water column in a good agreement with the experimental data and provide limiting condi- 
tions typical of the south part of Mediterranean Sea during the summer. Changes in the phytoplankton composition, 
both qualitatively and quantitatively, are related to the different depths considered, with light intensity and nutrient 
availability being the most important factors. Picophytoplankton demonstrated greater ability for photoacclimation 
than nano- and micro-phytoplankton**^^^^ . In fact, a higher contribution of picoeukaryotes to the phytoplankton 
biomass is ob seryed , specifically pelagophytes and prymnes iophyte s , which were also found to thrive elsewhere in 
cyclonic eddie^^^*^. This ability was also observed in culturd^SEIlIl, 

On the basis of our theoretical findings we can conclude that the position of the deep chlorophyll maximum depends 
on the parameter values used in the model. We used values of the buoyancy velocity v and vertical turbulent diffusivity 
-Df,, for which no oscillations occur. In this work we used the condition Df, = Dfi = 0.5 cm? /s, corresponding to poorly 
mixed waters along the whole water column, which causes the phytoplankton peak to have a width of few meters, as 
observed in the experimental data. Moreover, we also considered in our model the presence of an upper mixed layer, 
above the thermocline, characterized by a higher value of the diffusion coefhcients = Dn = 50 cm^/s), keeping 
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TABLE III: Results of x^, reduced chi -square [y^) goodness-of-fit test for site L1129b (left panel) and site L1105 (right panel) 
for different values of gr (stochastic dynamics - case 2). The number of samples along the water column is n = 176 for site 



L1129b and n = 563 for site L1105. 
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FIG. 12: Magnitude, depth, and width of the DCM as a function of an obtained from the model for site L1129b (panels a, b, 
c) and site L1105 (panels d, e, f). 



Db — Dfi — 0.5 en? I s for greater deptH^Sl. The results (here not shown) did not evidence any variations in the 
picophytoplankton distributions respect to the case of uniform diffusion coefficients = Dn = 0.5 err? j s) along 
the whole water column. This can be explained noting that in the ecosystem considered here the mixed layer, due to 
the depth of the thermocline, is not enough thick to influence the DCMs of the chlorophyll distributions. 
In our ecosystem the position and stability of the chlorophyll maximum, obtained from the model, depend not only on 
the vertical turbulent diffusivity, but also on the nutrient concentration at the bottom i?i„ and the maximum specific 
growth rate r. We also note that the values of Rin used in our model are c ompatib le with the nutrient concentrations 
measured along the water column in several sites of the Mediterranean SecPSESEll, 

Our numerical results were calculated by setting the maximum specific growth rate r at a value consistent with ex- 
perimental observations. Specifically, this value has been chosen so that the net per capita gro wth rat e g{z,t), used 
in the model, is in a good agreement with those experimentally observed for the picoeukaryote^sHzHini 

We recall that the estimations of the chl a content per picoeukaryote cell are highly variable, depending on the 
depth and water properties (oligotrophic or eutrophic) examined. Moreover these estimations reflect the taxonomic. 
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ecological and physiological diversity and the plasticity highlighted in previous studies I^^HISI^ jj^ our model we took 
into account this aspect. In particular, after obtaining the numerical results for the phytoplankton concentration 
expressed in number of cells/m^, we used the experimental findings given in Ref !^ to convert the numerical results 
into chl a concentration expressed in fJ.g/1. Specifically, because of the peculiarities of our model, suitable to describe 
the dynamics of the picoeucaryotes, we used the conversion curves typical of these species and compared the results 
with the experimental chl a concentrations sampled in two different sites of the Mediterranean Sea (Channel of Sicily). 
From the comparison we found that the values of chl a concentration obtained numerically are in a good agreement 
not only with our data but also with those measured by Brunet et al.^^. In addition, we note that our numerical 
results for the picoeu karyo te concentration expressed in number of cells/m'^ match the corresponding experimental 
data reported in Refs.BSHU 

More precisely, as a first step we used a deterministic model, consisting of an auxiliary equation for the light inten- 
sity and two differential equations, one for the dynamics of the phytoplankton biomass, the other for the dynamics 
of the nutrients. The numerical results showed a good qualitative agreement with the real data, even if discrepancies 
were observed between the characteristics of the chl a concentration profiles provided by the model and those obtained 
from the real data. 

To improve the agreement between numerical and experimental distributions, we modeled the random fluctuations 
of the environmental variables, by adding a term of multiplicative Gaussian noise in the differential equation for 
the phytoplankton biomass. The results obtained indicate that the presence of random fluctuations, acting directly 
on the phytoplankton biomass, determines chl a stationary distributions more similar to the experimental ones. In 
particular, we found that both the position and magnitude of the DCMs agree very well with the experimental find- 
ings. Afterwards, we modified the deterministic model considering the role of a noise source which influences directly 
the dynamics of the nutrients, by adding a term of multiplicative Gaussian noise in the differential equation for the 
nutrients. In this case we observed for suitable noise intensities (much lower than those used in the equation for the 
phytoplankton biomass) a further improvement of the numerical distributions of chl a concentration respect to the 
experimental ones. In addition, we found that higher noise intensities (comparable with those used in the equation for 
the phytoplankton biomass) , cause a rapid extinction of the phytoplankton community. The results obtained indicate 
that the proposed stochastic model is able to reproduce patterns of real phytoplankton distributions when aquatic 
ecosystems with poorly mixed waters are considered. 
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